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Abstract 

The static and dynamic critical properties of the ferromagnetic q-state Potts models on a square 
lattice with q = 2 and 3 are numerically studied via the nonequilibrium relaxation method. The 
relaxation behavior of both the order parameter and energy as well as that of the second moments 
are investigated, from which static and dynamic critical exponents can be obtained. We find 
that the static exponents thus obtained from the relaxation of the order parameter and energy 
together with the second moments of the order parameter exhibit a close agreement with the exact 
exponents, especially for the case of q = 2 (Ising) model, when care is taken in the choice of the 
initial states for the relaxation of the second moments. As for the case of q = 3, the estimates for 
the static exponents become less accurate but still exhibit reasonable agreement with the exactly 
known static exponents. The dynamic critical exponent for q = 2 (Ising) model is estimated from 
the relaxation of the second moments of the order parameter with mixed initial conditions to give 
z(q = 2) ~ 2.1668(19). 
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I. INTRODUCTION 

The so-called nonequilibrium relaxation (NER) method [l, 2] deals with the nonequilib- 
rium critical relaxation of a statistical model instantaneously brought to its critical temper- 
ature from its nonequilibrium initial state, typically a fully ordered state. A unique feature 
about the NER method is that it allows one to determine the equilibrium critical properties 
of statisitical systems from their nonequilibrium relaxation kinetics. Therefore, together with 
the method jj, 4] using a short-time critical dynamic properties f| (see Section III.D), 
the NER method has been an interesting and valuable tool to study the statistical systems 
whose critical properties are unknown. 

With recent advances in computing power, it seems worthwhile to reinvestigate in more 
depth the potentials or limitations of the NER method in terms of numerical accuracies in 
static and dynamic exponents. In this work we chose the ferromagnetic g-state Potts model 
|3] on a square lattice (with q = 2 and 3) as a testground for the NER method. It appears 
to us that in applying the NER method, the energy relaxation itself has been overlooked in 
evaluating the critical exponents: previous works a, |9| employed mainly the second moments 
involving the energy, ignoring the relaxation of the energy itself. In this work we utilized the 
critical energy relaxation as well as that of the order parameter for evaluating the critical 
exponents. Due to better self-averaging of these single moment quantities, higher accuracy 
is expected for the estimates of the exponents. Since for the ferromagnetic g-state Potts 
model the critical temperature and energy are exactly known for general values of q, and 
the static critical exponents are exactly known for q < 4 [3, Q , our work can be considered 
as a calibrational study on the NER method as to how accurate estimates for the critical 
exponents can be given via the NER method. 

We find that by combining the relaxation of the two first-moment quantities, i.e., order 
parameter and energy one can obtain the ratio of the equilibrium critical exponents, through 
which the nature of the phase transition or, equivalently, the universality class, may be 
determined regardless of our knowledge of the exact value of the dynamic exponent z. By 
incorporating the time dependence of the second moment of the order parameter (starting 
from either disordered or ordered initial states), the dynamic exponent can be independently 
obtained. This can be combined with the relaxation of the first moment quantities to 
give the two independent static exponents. We here find that the second moment of the 
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order parameter obtained with disordered initial state gives the static exponents with higher 
accuracy. Alternative way of obtaining static exponents is to employ the second moments 
involving energy fluctuations. However, the energy fluctuation appears to be influenced by 
strong logarithmic corrections in the case of q = 2, which makes it difficult to estimate the 
exponents with high accuracy. 

The present work is motivated by our interest in applying the NER method to the statis- 
tical systems whose critical properties are not completely understood. We are particularly 
interested in clarifying the long-standing controversy on the critical properties of the FFXY 



model Uj on a square lattic e [121] and the antiferromagnetic XY model 13] (or the antifer- 



romagnetic clock model 







., uM) on a triangular lattice. Even though there exist previous 
studies 0, [l6, 0, [3] from nonequilibrium dynamics perspective on this problem, we hope 
to undertake more careful studies on the FFXY models by applying the NER method with 
the strategy presented in this work. 



II. MODEL AND SIMULATION METHOD 



In the ferromagnetic g-state Potts model on a square lattice each spin <j{ can take one of 
the q possible values, <7j = 1, 2, • • • , q, and the spin interaction is defined by the following 
Hamiltonian 

H = -J J2 <K (J i> "i)> cri = l,2,---,g (1) 

<ij> 

where 5(a,b) is the Kronecker delta function and < ij > denotes the nearest neighbor spin 
pair. J(> 0) is the interaction strength. 

The model is known to exhibit a continuous phase transition for q < 4 and discontinuous 
ones for q > 4 in two dimensions For q = 2, the Potts model is equivalent with the Ising 
model since S(ai, <x,-) = (1 + SiSj)/2 with Si = ±1. It is remarkable that though the model is 
not exactly solvable except for q = 2, the critical temperature T c and energy E c are exactly 
known for general q in some two dimensional lattices. For the square lattice they are given 

byfl 

where is the Boltzmann constant. Here and after we set ks = 1 and J = 1. Hence 
T c = 1.13459 E c = -1.70711- •• for q = 2, and T c = 0.99497 E c = -1.57735- •• for 
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q = 3. The exact static critical exponents are given in TABLE I. 

We will be interested in characterizing the time evolution of the system instantaneously 
'heated' to T c from a fully ordered initial state. The system then proceeds toward the 
equilibrium via Monte Carlo kinetics with the Metropolis algorithm for the flip of a ran- 
domly selected spin. Since we are interested in the long time dynamics of the system in 
the thermodynamic limit, we take the system size large enough so that the simulation can 
appropriately mimic the nonequilibrium relaxation of the infinitely large system within the 
simulation time window. This is similar to the case of phase-ordering kinetics [3] in which 
the equilibration time is practically infinite due to large system size, making the time evolu- 
tion ever nonequilibrium. We use the square lattice with linear size L = 600. The maximum 
simulation time is t max = 10 4 Monte Carlo steps (mcs). The presented results are averages 
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over 20, 000 samples. Compared to the existing dynamic studies (4, 
the present model, our work has employed bigger lattice size with one or two decades longer 
Monte Carlo steps of simulations. Within the simulation time window, we have checked that 
there is no finite size effect on our simulation results. 



III. MEASUREMENTS AND DISCUSSIONS 

We first define the local order parameter rrii(t) and local energy ei(t) for the model, and 
their respective sums A4(t) and S(t) as 

1 1 (i) 

m ^ = (VT TyM^)>^(0)) - l), e,(t) = --E^iW.^iW) 

M(t) ee = (3) 

i i 

where Y^j denotes the sum over the nearest sites of the site i. We have set 0^(0) = 1 in most 
of the present simulations. In the NER method, the determination of the critical exponents 
involves the cumulants of A4(t) and S(t). The order parameter M(t) and the energy density 
E(t) are given by M(t) = (M{t))/N and E(t) = (£(t))/N respectively where N is the total 
number of spins, and the angular bracket < • • • > denotes average over different realizations 
of the time evolution. 



4 



A. Locating the critical temperature and energy 

The NER method provides us with an efficient tool to determine T c with high accuracy. 
If the temperature T is above or below T c , then the order parameter would relax (stretched) 
exponentially in time to its equilibrium value which is zero for T > T c , or nonzero for 
T < T c in the continuous transition as in the present model. The relaxation time scale r 
depends on how close T is to the critical temperature T c . Approaching T c , r exhibits a 
power law divergence. Below T c , the order parameter M(t) relaxes toward the nonvanishing 
equilibrium value. Since only at T c the order-parameter exhibits a power-law relaxation, 
T c can be located as the temperature at which the order-parameter relaxation gives the 
best straight line in a log-log plot of M(t) versus t. More accurate location of T c naturally 
requires longer simulation time. 

Since T c is exactly known for the present model, we can provide a test example how the 
order-parameter relaxation can be used to narrow down the critical temperature, as shown 
in Fig. 1(a) (q = 2) and Fig. 1(b) (q = 3). Figure la shows the relaxation of the Ising 
(q = 2) order parameter for various temperatures near T c . M(t) shows strong upward and 
downward curvatures at T = 1.13 and at T = 1.14, respectively. So the critical temperature 
should be in the range 1.13 < T c < 1.14. We find in this way that the temperature range can 
be readily narrowed down to ST = 4 x 10~ 4 within the present simulation time (t max = 10 4 
mcs). That is, we have 1.1342 < T c < 1.1350. Higher accuracy would be achieved with 
longer simulation time. In the same manner, for q = 3, T c can be easily located within the 
range ST = 5 x 1CT 4 as demonstrated in Fig. 1(b). 

In contrast to the case of the present system, the exact value of the critical energy is 
not known for many other systems. For such cases, using an ansatz of the critical energy 
relaxation at T c , one can determine E c by tuning the value of E* such that (E* — E(t)) 
versus t gives the best straight line at long times in the log-log plot. Figure 2 shows such an 
example for the present system. In our simulation time window, one could readily narrow 
down E c within the range SE ~ 10~ 3 . 
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B. Relaxations of the first moments 



We first measure the relaxation of the order parameter M(t) and the energy density 
E(t). When the system is instantaneously brought to T c starting from a fully ordered state, 
both quantities exhibit power-law relaxations toward their equilibrium values si, 0, 
For better analysis of the relaxation behavior and estimation of the critical exponents, we 
included correction terms to the leading-order scaling behavior in the following form (see 
APPENDIX for a derivation of the leading-order critical relaxation) 

M(t) ~ t-^(l+^ + ^L + ^L+. • •), E c -E(t) ~ t-^-^(l + ^+^+$+. ■■) 

V ' \ t &M t 5 M t & M > V ' V t &E t 5 E t 5 E J 

(4) 

where z is the dynamic exponent, /3 and v the static exponents, and d the spatial dimension 
{d = 2 in the present case). In (4), cm, c' m , c" m and ce, c' e , c" e are constants, and 5m, $e, etc 
are the exponents of the correction-to-scaling terms. Simulation results for the relaxations 
of M(t) and (E c — E(t)) are shown respectively in Fig. 3(a) and Fig. 3(b). 

It is easy to see that the general form of Eq. (4) with all the correction exponents kept 
independent, is not appropriate for fitting. For example, if we put all the correction expo- 
nents equal to one another, then we can generate the same function with many different 
combinations of the coefficients provided that the sum of the coefficients are the same. This 
will cause the fitting procedure with a general initial guess to fail to converge to unique 
stable values for the values of exponents and other coefficients. Therefore, we should impose 
some constraints on the correction exponents to achieve a stable fit. One of the schemes 
we tried is to let the first correction exponent to be a free fitting parameter but the n-th 
correction exponent are constrained to be equal to the first correction exponent plus (n — 1). 
That is, we put 5' M = 5m + 1 and 5" M = 5m + 2, etc (which we call FITnA, where n refers 
to the number of the correction terms). We also tried a scheme where the first correction 
exponent to be a free fitting parameter but the n-th correction exponents are constrained to 
be equal to n times the first correction exponent. That is, we put 5' M = 25m and 5% = 35m, 
etc (which we call FITnC). A special scheme of correction is to assume an analytic form 
for the correction terms where the correction exponents are all constrained to integer values 
(which we call FITnB). The cases of 5m = 1 in FITnA and FITnC reduce to the scheme of 
FITnB. The special case of fitting with a simple power law without correction terms will be 
referred to as FITO. 
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First consider the case of q = 2. In this case, we incorporated up to third correction terms 
(i.e., FIT3A, FIT3C and FIT3B) in our fit. In the cases of FIT3A and FIT3C, for reasonable 
convergence of the fitting procedure and estimation of the errorbars, we tuned the value of 
the coefficients of the highest order correction terms (i.e., d' M and c" E ). These values of c" M 
and c" E were tuned (chosen) such that the resulting fit function shows good agreement with 
the data. For some finite range of these coefficient values (c' M ranging from 0.0 to 0.1), the 
quality of the fit is found to be acceptible. The median of the values of the exponents was 
chosen as the representative value of the exponent and the range of the values can be used 
to estimate the error bars. As for the relaxation of the magnetization, we found that, as the 
value of c" M was varied, the value of the dominant relaxation exponent (3/zv showed only 
small variation. We also found that, in terms of the value of the dominant exponent (3 / zv 
and the first correction exponents, the two schemes of FIT3A and FIT3C gave almost the 
same results (see TABLE II). If we take the median of the fitted values from this analysis 
as the most probable value of (3/zv, we get 



Now the corresponding value of the first correction exponent 5m from the above fit exhibits 
a monotonic increase from 0.91 to 1.06 (FIT3C) and also from 0.90 to 1.07 (FIT3A) as the 
value of c' M increases from 0.0 to 0.1 (TABLE II). Taking the median of the combined range 
(from 0.90 to 1.07) we may take the representive value of 5m as 5m — 0.985(85). This may 
indicate that the exact value of the first correction exponent is equal to one. Incidentally 
we also tried to fit the data with analytic corrections in which the first correction exponent 
is set to be equal to one, the second correction exponent two, etc. We found that the fitting 
quality was excellent with essentially the same dominant exponent as the representative 
value given above. 

As for the relaxation of the excess energy, we also applied the above fitting schemes. 
However, in this case of excess energy, we had to directly tune the value of the first correction 
exponent in order to obtain a reasonable fit to the data. We found that the best fit could 
be obtained in the range of the value of 5e between 0.90 and 1.05 (see TABLE III) which 
again is in reasonable agreement with the case of the magnetization relaxation (Fig. 3(b)). 
The values of the fitted relaxation exponents are 



0.057722(42) 



for q = 2. 



(5) 



vd- 1 



0.46227(71) 



for q = 2. 



(6) 



zv 
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Now we turn to the relaxation dynamics of q = 3 model. As for the magnetization 
relaxation, we tuned the value of the coefficient c' M or d' E of the highest order correction 
terms. We find that both schemes FIT3A and FIT3C give approximately the same results 
for values of the relaxation exponent (3/ zv (TABLE IV). Combining the results of FIT3A 
and FIT3C we obtain the values of the exponent as 

— = 0.06020(17) for q = 3. (7) 

zv 

Correspondingly the value of the first correction exponent 5m ranges between 0.65 and 0.98. 
As for the relaxation of the excess energy, we also applied the above fitting schemes. We 
found that the best fit could be obtained in the range of the value of Se between 0.6 and 
0.74 which exhibit a narrower range for the correction exponent compared with the case of 
magnetization relaxation. The value of the fitted relaxation exponent for excess energy is 
(TABLE V) 

vd ~ l = 0.3626(40) for q = 3. (8) 

zv 

Figures 3(c) and 3(d) show the time dependence of M A (t) = M(t)A t l3 ^ zu for q = 2 and q = 3 
respectively together with the fitting functions. Here Aq l t~^' zv represents the asymptotic 
scaling obtained from the fitting. Also, Figs. 3(e) and 3(f) show the time dependence of 
E A (t) = (E c - E(t))B & d - 1 ^ zv with the fitted values of the exponents for q = 2 and q = 3 
respectively together with the fitting functions, demonstrating the asymptotic nature of the 
scaling behaviors. Here BQ 1 t~v y ' d ~ l >' Z1 ' represents the asymptotic scaling obtained from the 
fitting. 

Note that the value of 13/ zv is particularly small. This implies that the order-parameter 
relaxes much more slowly than the energy does: even after four decades of relaxation time, 
the order parameter relaxed only half of its initial value. This is another reason why the 
statistics is better for the order parameter relaxation than for the energy relaxation. The 
above values of (3 / zv for both q = 2 and 3 can be compared with the results obtained from 
the short-time dynamics both on square % |2l|] and triangular 25j lattices. 

Since (3/ zv and (yd — l)/zv contain the common factor 1/zv which involves the dynamic 
critical exponent z, taking the ratio of the former to the latter eliminates the exponent z, 
yielding (3/(vd — l) which involves only the two static exponents. Then the numerical values 
of the ratio (3/(vd — 1) obtained from the above fit can be compared with the corresponding 
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exact values as 



0.1249(4), 1=0.125 for q = 2, 



vd - 1 W - 1/ exact 

0.1660(23), (—r—r) f = - = 0.16667 for g = 3. (9) 



i/d- 1 "" W- 1/ exact 6 

We see that for both cases of q = 2 and q = 3, the values of (3/{vd — 1) are quite close to 
the exact theoretical values. 

Elimination of the dynamic exponent z in the above ratio is an important feature that 
may be used to distinguish the universality class of systems whose critical properties are not 



known. Prominent examples of such systems are the two-dimensional FFXY models [12 ]. 
We point out that as is done here, measuring the relaxations of the chirality order parameter 
and the energy for such systems one can determine whether the chirality transition in that 
system belongs to the Ising or the 3-state Potts universality class, or to another universality 
class, without knowing all the static exponents separately. 

One can obtain an estimate for z from the fitted exponents for the relaxation of mag- 
netization and energy as given above by making use of the fact that the static exponents 
are exactly known in the case of q = 2 and q = 3. This is perhaps one of the easiest ways 



of obtaining z for systems whose static critical properties are exactly known [28l. l29j|. If we 
use the magnetization relaxation, the values for z derived from the values of the relaxation 
exponents in (5)-(8) (assuming the exact theoretical values of the static exponents) are given 
by 

z M = 2.1656(16) z E = 2.1632(33) for q = 2, 

z M = 2.2150(62) z E = 2.2064(110) for q = 3, (10) 

Here values of the exponents are the median values in the optimal fitting regime, while the 
error bars are estimated from the maximal dispersion of the fitted values. 



C. Relaxations of the second moments 



The above procedures of obtaining z of course cannot be carried out for systems with 
unknown critical properties. It is thus desirable to provide ways of independently measuring 
the dynamic exponent. We present below one systematic way of doing it by considering the 
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second moments o 
behavior in time 



the order parameter, C MM.it) i which exhibits a leading-order scaling 
as (see APPENDIX) 

C MM (t) = N\(M\t)) - (M(t)) 2 } ~ (11) 



Note that (11) is valid for disordered initial states as well. The dynamic exponent z can be 
isolated if one considers a time-dependent second moment /mm^) defined as 

lMM(t) = ~ ** (12) 

which is obtained using the leading-order behavior for M(t) given in (4). One usually 
measures fMM.if) i n two ways: one measures fMM{t) for fully ordered initial states, or one 
can first measure Cmm{^) f° r disordered initial states, and then divide it by M 2 (t) measured 
for the fully ordered initial states. The latter way is therefore using the two different types 
of initial conditions. We denote the former result by Jmm orderedit), and the latter one by 

n 

fMM,mixed(t). Alternative way [4j of obtaining z using the short-time critical dynamics is 
described in the next section. 

Shown in Fig. 4(a) and Fig. 4(b) are fMMif) for g = 2 and q = 3, respectively. Here, 
we apply similar methods of analysis as in the previous section. That is, fMM(t) ~ 

t d/z A + CMM Jt s MM _|_ y We found that for the case q = 2 (using both FIT2A and 

FIT2C) f M M,ordered{t) and f M M,mixed(t) give (TABLE VI and VII) 

Zordered = 2.1545(52), z mixed = 2.1668(19) for q = 2. (13) 

We see that there exists small but non-negligible discrepancy between the exponents z mixec [ 
and z ori iered- We also note in particular that z m i xe d is very close to zm obtained from the 
relaxation of the magnetization in (5) (assuming the exact values of the static exponents). 
Figure 4(c) shows the time dependence of the quantity fMM,A{t) = fMM{t)Cot~ d l z for the 
case of mixed initial states (for q — 2) together with the fitting function, where C^ x t d l z 
represents the asymptotic scaling obtained from one of the best fits. 
In the case of q = 3, z or d er ed and z mixe d exhibit a larger discrepancy: 

^ordered = 2.1334(35), z mixed = 2.1735(40) for q = 3 (14) 



As first pointed out by Zheng J4], this difference between the two estimates is not due to 



the statistical error of the data. This difference is also observed in triangular lattice 
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25] 



More severe disagreement between z or dered and z mixe d was reported for the Baxter- Wu model 



32 



33] : z ordered = 2.07(1) [34J and z mixe d = 2.294(6) 



351 ] . We suspect that the underlying 



reason for this discrepancy is that the time scale at which the scaling behavior sets in, may be 
much longer for fMM,ordered(t) than that for fMM,mixed{t) as q gets larger, which may be due 
to the nature of the broken symmetry of the initial states related to the higher degeneracy 
of the ground states for larger q. It is thus expected that the stronger discrepancy will be 
observed for the 4-state Potts model as well. 



1. The static exponents v and (3 

The values of dynamic critical exponents obtained from second moments of the order 
parameter can now be combined with the relations of exponents obtained from the relaxation 
of the first moments of the order parameter and energy to give the static exponents such as v 
and (3. Since it appears that the value of z m i xe d (especially for q = 2) is more consistent with 
the results of magnetization and energy relaxation, we substituted z mixe d for the dynamic 
exponents in the relaxation of magnetization and energy. For example in the case of q = 2, 
substituting the value of z mixe d = 2.1668(19) in (13) for z in the relations (5) and (6) gives 
\jv — 0.9984(25) and (3 = 0.12527(51). These values yield close estimates to the exact 
exponents: \/u = 1 and (3 — 1/8 = 0.125. 

As for q = 3, when the value of z m i xe d = 2.1735(40) in (14) is substituted for z in the 
relations (7) and (8) we obtain \ jv = 1.212(16) and (3 = 0.1080(20), which appears to exhibit 
approximate agreement with the exact values \ jv = 6/5 = 1.2 and (3 = 1/9 = 0.1111. But 
this is definitely less accurate than the case of q = 2. If we use z or dered instead of z m i xe d, then 
we obtain even less accurate results for the static exponents, which is naturally expected 
from the discrepancy of the value of z or dered from z obtained from the relaxation of order 
parameter or energy assuming the exact static exponents. One thus should be cautious in 
determining z from the second and higher moments of the order parameter for systems with 
unknown critical properties. 

We have so far followed the procedure of using the relaxations of the order parameter, 
the energy, and the second moment of the order-parameter. An alternative way js], [9} of 
obtaining the exponent v is to use the second-order moments involving energy fluctuations, 
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fMs(t) or fee(t), which exhibit the leading-order scaling behaviors as [8J (see APPENDIX) 



f M s(t) = N 
fse(t) = N 



< M(t)£(t) > - < M(t) >< £(t) > 

<M{t) ><£(t) > 

< £ 2 {t) > - < 8{t) > 



fl/zv 



t a/zu (15) 



< £(t) > 2 

For q = 2, we first obtain the leading-order exponent 1/zu = 0.4624(29) by using FIT3A 
(Fig. 5(a)). However, we are unable to extract the exponent (2 — vd)/zv from fss{t) since 
it exhibits a strong curvature with the local slope decreasing in time, as shown in Fig. 5(b). 
This may indicate that the exponent (2 — dp) eventually vanishes in the long time limit, 
implying v = 1. We suspect that this strong curvature in the log- log plot of fse(t) may 
correspond to a logarithmic correction. Figure 5(c) shows such a fit with fgs(t) ~ (logt)* 
with <p ~ 0.80. From the value of the above exponent 1/zu = 0.4624(29) we obtain 1/v — 
1.0019(72) using z mixe d = 2.1668(19). This estimate compares quite well with that obtained 
from M(t), E(t), and f M M(t). 

For q = 3, in contrast to the Ising (q = 2) case, both moments fMsif) and fse(t) exhibit 
asymptotic power-law behavior. The average slopes obtained from Fig. 5(a) and Fig. 5(b) 
are given by 1/zu = 0.5515(61) and (2 — ud)/zu = 0.1911(72). These values yield the static 
and dynamic exponents 1/v = 1.210(58) and z = 2.193(129). It is interesting to see that 
this dynamic exponent is larger than z m i xe d = 2.173(2) given in (14), and is rather closer to 
those given in (10). However, one should also note that the error bars are much larger in 
this case. 

One can therefore conclude that for the present system using the relaxation of the order 
parameter and the energy combined with the value of z m i xe d yields better estimates for static 
exponents than using the second-order moments involving the energy fluctuations. From this 
point of view, the moments fMeif) an d fss{t) may be used to check the consistency of the 
results obtained using the energy relaxation itself. 



D. The short-time critical dynamics 

At this point, for comparison, it is worthwhile to briefly describe the method using the 
short-time critical dynamics which offers an alternative way of determining the static and 
dynamic exponents. More details can be found in Q]. When quenched to T c from the 
initial disordered state but with sufficiently small 'magnetization' M , the system exhibits a 
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nonequilibrium scaling regime in which the order parameter typically shows an anomalous 
power-law increase as 

M(t)~M t 9 (16) 

before eventually relaxing toward its equilibrium value, where the exponent 9 is a new 
nonequilibrium critical exponent. Though it typically takes a positive value, it becomes 
negative for some systems such as the tricritical system [38(, the Ashkin- Teller model 
the Blume-Capel model 40] , the 4-state Potts model 41[ , and for the BW model 



34 



35 



39 
42 



The time scale tm associated with the power-law increase depends upon the magnitude of 
the initial magnetization as tm ~ Mq Z ^ x ° where the scaling dimension xq is related to the 
other exponents clS Xq — 9z + 0u. 

In practice, 9 was obtained by linearly extrapolating the effective exponents 9 e ff (for 
finite Mo) vs. M to M = 0. Alternative way of obtaining 9 was first proposed by Huse Q]. 
Using a scaling ansatz for the nonequilibrium structure factor of the two-dimensional kinetic 
Ising model, Huse has shown that the auto-correlation function of the global order-parameter 
M. (t) with random initial conditions exhibits a power-law behavior 



A M (t) = - < M(t)M(0) >~t e 



(17) 



Huse has also shown that the spin auto-correlation function A(t) exhibits a critical decay in 
time as 



A(t) 



q 1 
q - IN 



< E (*(*(°W*)) 



>~ t~ 



r 



(18) 



Tome and Oliveira [43J] have proved that (18) holds for the equilibrium as well as nonequi- 
librium systems (he. without detailed balance) possessing the global up-down symmetry. 
Later on, Tome [44j has shown that the validity of (18) is extended to systems with other 
symmetries. It is much more convenient to measure 9 by measuring the auto-correlation 
Am{^) since one is free from cumbersome preparation of the initial states with very small 
magnetization and an extrapolation procedure to vanishing initial magnetization. 

The new nonequilibrium exponent 9 has been obtained as follows. For q = 2, 9 = 0.191(1) 



(Heat Bath) and 9 = 0.197(1) (Metropolis) using (16) |4|, [23|, [2j, [25| , and 9 = 0.19 using 
(17) [6]. In addition, Grassberger [45[] also measured the exponent 9 as 9 = 0.191(3) using a 



damage spreading method. Likewise, for q = 3, The exponent 9 takes t 
9 = 0.075(3) (Heat Bath) and 9 = 0.070(2) (Metropolis) using (16) 



re following yalues: 



23|, 



24 



25|, and 
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9 = 0.072(1) using (17) [4l|. The exponent 9 was also measured using (17) as 9 = 0.093(4) 



46j for a nonequilibrium cellular automaton model with C$ v symmetry which belongs to 



the 3-state Potts universality class 47] . In addition, there exists measured value of 9 for 
q = 4 using both (16) and (17), which are given respectively by 9 = —0.0471(33) and 9 = 
—0.0429(11) 41]. It is an interesting fact that the 4-state Potts model and the BW model, 
which belong to the same static universality class, do not share the same nonequilibrium 
exponent 9: 9 = -0.185(2) [42] or -0.186(2) [35J] for the BW model. It would be interesting 
as well to clarify the question as to whether this nonuniversal dynamic behavior is extended 
to the case of z in both systems. 

Independent measurements of the exponents 9 and A using (17) and (18) provide an 



alternative way 



20 



d/(X + 9). The measured 
2.155(3) for q = 2, and 



23j of obtaining z through the relation z 
values of A and hence z are given by A = 0.737(1) and z = 
A = 0.836(2) and z = 2.196(8) for q = 3 [4|. As for the BW model, A = 1.188(10) and 
z = 1.994(24) 42j were obtained using the continuous-time Monte Carlo algorithm. One 



disadvantage of this method is that the rather large value of A induces strong fluctuations 
in its measurement, which is the main source of the error for the estimate of z. 

With z in hand, one can use (11) with disordered initial states to obtain the ratio /3/v, or 
equivalently, the exponent r] = 20/ v (for d = 2). In order to determine (3 and v separately, 
one can use the following scaling ansatz for the order parameter near and below T c 



M(t, e) = r^(^ w e) 



(19) 



where JF is the scaling function, and e is the reduced temperature defined as e = (T c — T)/T c 
The equation (19) reduces to (4) for e = 0. Then one can give an estimate for 1/v using 



d_ 
d~e 



lnM(t,e)| e=0 = t 1/ ^(^ln^(x 



e=0 



(20) 



jv — 1.027(6) (triangular lattice) 



This method gives 1/v = 1.03(2) (square lattice) [l] and 

2M for o = 2, and 1/v = 1.24(3) (square lattice) |j, bl| and 1/v = 1.223(8) (triangular 
lattice) [25] for q = 3. But, involving the difference of the order parameter at several close 



temperatures near T c may cause considerable error in the estimate of the exponent v. Our 
method (pesented in Section III.C) utilizing the energy relaxation can provide an alternative 
and better way of estimating the exponent v without using (20). 
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IV. SUMMARY AND CONCLUDING REMARKS 



Through Monte Carlo simulations we have investigated the nonequilibrium critical dy- 
namics of the ferromagnetic g-state Potts model on a square lattice. Primary purpose of the 
present work was to evaluate how accurate estimates the NER method can provide for both 
static and dynamic exponents of the statistical systems. The ferromagnetic Potts model is 
an ideal test system for that purpose since the static critical properties of the model are 
exactly known. 

In contrast to other nonequilibrium approaches, we utilized both the order parameter 
and energy relaxation. The ratio of the power-law exponents of the order parameter and 
the energy involves solely the static exponents. In order to separately measure the two 
independent static exponents, one first measures the dynamic exponent z by considering the 
order-parameter moments starting from either random or fully ordered initial states. It is 
found that in the present model, when the fully ordered initial states are employed, the time 
scale associated with the asymptotic scaling for the relaxation of the second moments of the 
order parameter sets in, seems to be much longer for higher q than that for the disordered 
initial states and may exceed the present simulation time window. As a result, the slopes 
obtained from the second moments of the order parameter may differ. One thus has to be 
very careful in obtaining the dynamic exponent z from the time-dependent moments for the 
order parameter. Once z is determined, the two independent exponents v and (3 can be 
determined from the relaxation of the energy and the order parameter, respectively. 

The present work has demonstrated that, in the case of q = 2 (Ising) model this method 
can provide accurate estimates for the static and dynamic exponents provided that, as 
mentioned above, special care is taken in the choice of the initial states for the estimation 
of the dynamic exponent z. For example the dynamic exponent z m i xe d obtained from the 
relaxation of the second moment of the order parameter with disordered initial states are 
more consistent with the relaxation of the order parameter and excess energy when the 
exact values of the static exponents are assumed. In the case of q = 3 Potts model, even 
though the static exponents obtained with the same methods are less accurate than in the 
case of q = 2, they still show reasonable agreement with the exact known values of the static 
exponents. 

The fact that the dynamic exponent z does not appear in the ratio of the power-law 
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exponents of the order parameter and energy, may be used to clarify the nature of the phase 
transition whose critical properties are under debate. An outstanding example is the phase 
transition associated with the chirality order in the FFXY models. The relaxations of the 
order parameter and energy are expected to become slower for frustrated systems. Hence 
the above ratio can be measured with higher accuracy. We therefore tend to believe that 
once the critical temperature is accurately determined, the critical relaxations of the order 
parameter and energy for the chirality phase transition can tell whether the phase transition 
belongs to the Ising universality class, or to the 3-state Potts class, or else to an entirely 
different universality class. 
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APPENDIX: CRITICAL RELAXATIONS OF MOMENTS OF THE ORDER PA- 
RAMETER AND ENERGY 

The critical relaxations of the order parameter and energy, and their higher moments can 
be obtained from the following scaling behavior of the nonequilibrium generating function 
ip(e, h,t) which reduces to the equilibrium free energy density in the limit t — > oo: 

tfj(e,h,t) = b- d ^(e-W T ,h-b y \t-b- z ) 

= r d/z ^(e ■ t yT/z , h ■ t yh/z , 1) (A.l) 

where b is the scaling factor, e = \T C — T\/T c the reduced temperature difference, and h the 
external (magnetic) field coupled to the order parameter. The above result was first derived 



by Suzuki 48]. In (A.l), y? and y^ are the scaling dimensions associated respectively with 
temperature and external field, which are given by yx = 1/V and yh = d — (3/u. The 
relaxation of the order parameter M(t) at T c is obtained by differentiating the generating 
function ip(e, h, t) with respect to the magnetic field h: 

M(t) = (^) e _ h _ Q ~ t {yh - d)/z = r p f' v (A.2) 
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Likewise, the energy relaxation is obtained from the derivative of the generating function 
with respect to the reduced temperature e: 

(E(t) - E c ) = (^) _ h _ o ~ t(w-'0/* = r (»d-i)/uz = t -(i-a)/ Z u (A.3) 

where the hyperscaling relation 2 — a = ud was used in the last equality. The equations 
(A. 2) and (A.3) gives the leading order part of (4). 

The second order moment of the order-parameter is accordingly obtained by differenti- 
ating twice ifj(e, h, t) with respect to h: 

< (SM{t)) 2 >=< M 2 {t) > -M\t) = (^) e=h=Q ~ t^-^ z = t^ 2 ^ z ( A .4) 

The equation (A.4) gives (11). 

Likewise, the second-order moment of the energy is given by 

< (5£(t)) 2 >= (^) e _ h _ Q ~ t {2yT ~ d)/z = f( 2 -"«0A" = t a/zu (A.5) 

The equation (A.5) gives the last member of (15). 

Finally, the cross-moment of the order parameter and energy is given by 

< 6M(t)6£(t) >= (|^) e=/i=0 ~ t {VT+Vh - d)/z = t^ )/zu (A.6) 
Dividing < 5M(t)5S(t) > by M(t) gives 



< 6M(t)6S(t)> 1/zv 
M(t) 

The equation (A. 7) yields the first member of (15). 
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FIGURE CAPTIONS 



Fig. 1. The relaxations of the order parameter M(t) for (a) q — 2 and for (b) q — 3 at 
temperatures near T c ((a) T = 1.130, 1.1323, 1.1342, 1.13459(T C ), 1.135, 1.137, 1.14, 
and (b) T = 0.993, 0.994, 0.9945, 0.99497(T C ), 0.9955, 0.996, 0.997). The critical 
temperature can be determined as the temperature at which M(t) exhibits a critical 
relaxation (represented by a broken line) through a change from an upward curvature 
to a downward curvature. 

Fig. 2. The relaxations of the energy difference (E* — E(t)) in a log- log plot for (a) q — 2 
and (b) q = 3 with various test values of E* ((a) E* = -1.704, -1.705, -1.706, 
-1.70711(£ c ), -1.708, -1.709, -1.710, and (b) E* = -1.571, -1.573, -1.575, 
- 1. 57735 (E c ), -1.579, -1.581, -1.583). By tuning E* such that (E* - E(t)) ex- 
hibits the best power-law relaxation, one can determine both E c and — (ud — l)/zu. 
For each figure, the broken line represents the relaxation with the exact critical energy 
E c . 

Fig. 3. (a) The critical relaxation of the order parameter M(t), and (b) the critical relaxation 
of the energy difference (E c — E(t)). In (a) and (b) dot-dashed lines are curves obtained 
from FIT3A with optimal fits. Shown in (c) and (d) are plots of M A (t) = M{t)A^l zv 
(filled square) together with the fitting function (solid line) for q = 2 and q = 3 
respectively. Here Aq x t~^l zv represents the asymptotic scaling obtained from the fit- 
ting. Shown in (e) and (f) are the plots oiE A {t) = (E c - E(t))B & d -^/ zl/ vs. t (filled 
square) together with the fitting function (solid line) for q = 2 and q = 3 respectively. 
Here f?" 1 ^^/^ 1 )/ 21 ' represents the asymptotic scaling obtained from the fitting. 

Fig. 4. The order-parameter cumulant f'MM{t) vs. t in a log-log plot for q — 2 (a) and q = 3 
(b). (c): /mmA*) = fMM(t)C t- d / z vs. t in the case of mixed initial states for q = 2 
(filled square) together with the fitting function (solid line). Here Cq 1 ^ 2 represents 
the asymptotic scaling obtained from the fitting. Dotted lines are the straight lines 
indicating the asymptotic power law behavior. 

Fig. 5. The second-order moments fMs(t) ( a ) an d fee CO (b) vs. t in a log-log plot for q = 2 
and q = 3. (c) fee(t) vs. In t in a log-log plot for q = 2. This plot shows a logarithmic 
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time dependence of fee(t) ~ (log£)^ with ~ 0.80 for the long-time region (t > 100 
(mcs)). 
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z 


l/u (3 


77 = 2(3 jv (d = 2) 


a = 2 — vd 


q = 2 


present work 


2.1668(19) (zmm) 


0.9984(25) 0.12527(51) 


0.2501(16) 


-0.0032(50) 




short-time dynamics 


2.155(3) 


1.03(2) 0.1165 


0.240(15) 


0.058(38) 




exact 




1 1/8=0.125 


1/4=0.25 







present work 


2.1735(40) {ZMM, disordered) 


1.213(6) 0.1080(20) 


0.2618(83) 


0.350(22) 


q = 3 


short-time dynamics 


2.196(8) 


1.24(3) 0.1085 


0.269(7) 


0.387(40) 




exact 




6/5=1.2 1/9=0.1111 


4/15=0.2666 


1/3=0.3333 



TABLE I: The measured dynamic and static critical exponents. 
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r" 

JV1 


P/vz 


8m 


zm 




0.0 


0.0576815 


0.8956 


2.16707 




0.01 


0.0576911 


0.9165 


2.16671 




0.02 


0.0577003 


0.9366 


2.16637 


FIT3A 


0.03 


0.0577094 


0.9559 


2.16603 




0.04 


0.0577181 


0.9745 


2.16570 




0.05 


0.0577266 


0.9924 


2.16538 




0.07 


0.0577427 


1.0265 


2.16478 




0.1 


0.0577652 


1.0738 


2.16393 




0.0 


0.05768 


0.907 


2.1671 




0.01 


0.057689 


0.9264 


2.16679 




0.02 


0.057699 


0.9442 


2.16642 




0.03 


0.057708 


0.9612 


2.16606 


FIT3C 


0.04 


0.057717 


0.9774 


2.16573 




0.06 


0.0577338 


1.0078 


2.16511 




0.08 


0.0577494 


1.0363 


2.16453 




0.1 


0.0577642 


1.0633 


2.16397 


FIT3B 


0.07834 


0.05769 


1.0 (fixed) 


2.16675 



TABLE II: The dominant exponents and the first correction exponents for the relaxation of mag- 
netization, obtained from the two fitting procedures FIT3A and FIT3C for the case of q = 2. zm 
denotes the value of the dynamic exponent derived from the value of the second column for (3/zv 
using the exact values of the static exponents for (3 and v. c M represents the chosen values of the 
coefficients of the highest order correction terms for fits. For comparison, we also added (the last 
line) the result of a fit with FIT3B where the first correction exponent 8m is set equal to unity 
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Se 


(vd — l)/vz 


E 


ze 




0.90 


0.46201 


0.2028 


2.1645 


FIT3A 


0.95 


0.46227 


0.47124 


2.1632 




1.0 


0.46253 


0.76616 


2.1620 




0.93 


0.462169 


0.26146 


2.1637 




0.95 


0.462265 


0.40190 


2.1633 


FIT3C 


0.96 


0.462315 


0.47326 


2.1630 




0.97 


0.462367 


0.54537 


2.1628 




1.0 


0.462530 


0.766162 


2.1620 



TABLE III: The dominant exponents and the first correction exponents for the energy relaxation, 
obtained from the two fitting procedures FIT3A and FIT3C for the case of q = 2. Note that 
ze denotes the value of the dynamic exponent derived from the value of the second column for 
{vd — l)/vz using the exact values of the static exponents for v. Here the exponents 5^'s represent 
the chosen values of the first correction exponents for fitting. 
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r" 

M 


Blvz 


8m 


zm 




0.03 


0.0600345 


0.6415 


2.2209 




0.06 


0.0601009 


0.7063 


2.2185 




0.07 


0.0601192 


0.7260 


2.2178 


FIT3A 


0.08 


0.0601374 


0.7451 


2.2171 




0.09 


0.060154 


0.7634 


2.2165 




0.10 


0.060170 


0.7811 


2.2159 




0.12 


0.0601995 


0.8149 


2.2149 




0.15 


0.060241 


0.8630 


2.2133 




0.06 


0.060141 


0.7713 


2.2170 




0.09 


0.060183 


0.8139 


2.2155 




0.10 


0.0601959 


0.8274 


2.2150 


FIT3C 


0.12 


0.0602175 


0.8524 


2.2142 




0.15 


0.060249 


0.88845 


2.2130 


FIT3B 


0.244 


0.0603559 


1.0 (fixed) 


2.2091 


FITO 




0.06058 




2.2010 



TABLE IV: The dominant exponents and the first correction exponents for the relaxation of mag- 
netization, obtained from the two fitting procedures FIT3A and FIT3C for the case of q = 3. c M 
represents the chosen values of the coefficients of the highest order correction terms for fits. Also 
shown in the last two lines are the results of fits using FIT3B and FITO (with no correction terms) 
respectively. 
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Se 


(vd — 1) Ivz 


r" 

O 771 

£/ 


Z E 




0.55 


0.36008 


0.3898 


2.2218 


FIT3A 


0.60 


0.36142 


0.6459 


2.2135 




0.61 


0.36167 


0.6997 


2.2119 




0.62 


0.36192 


0.7544 


2.21045 




0.63 


0.36216 


0.8099 


2.20899 




0.64 


0.36239 


0.8664 


2.20756 




0.70 


0.36298 


0.7896 


2.20397 


FIT3C 


0.72 


0.36336 


0.9676 


2.20165 




0.73 


0.36356 


1.057 


2.2005 




0.74 


0.36375 


1.1464 


2.1993 



TABLE V: The dominant exponents and the first correction exponents for the energy relaxation, 
obtained from the two fitting procedures FIT3A and FIT3C for the case of q = 3. Here the 
exponents 8e represent the chosen values of the first correction exponents for fits. 





$MM 


r" 


d/z 


ZMM 




1.35 


0.0 


0.92323 


2.16631 




1.3 


0.0 


0.92311 


2.16659 




1.25 


0.0 


0.92298 


2.16689 


FIT3A 


1.20 


0.0 


0.92285 


2.16720 




1.15 


0.0 


0.92271 


2.16753 




1.0 


0.6923 


0.92314 


2.16652 




1.0 


0.0 


0.92224 


2.16863 


FIT3C 


1.34595 


0.0 


0.92335 


2.16603 



TABLE VI: The dominant exponents and the first correction exponents for the relaxation of the 
second moment of magetization with random initial states, obtained from the two fitting procedures 
FIT3A and FIT3C for the case of q = 2. Here, the fitting is mostly done up to second corrections. 
c!' M is the chosen values of the coefficients of the highest order correction terms for fits. 
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5mm 


Ivl lVl 


d/z 


zmm 

1V1 1V1 




0.8317 


0.15 


0.927756 


2.15574 




0.8471 


0.09 


0.927876 


2.15546 




0.8545 


0.06 


0.927933 


2.15533 


FIT2A 


0.8618 


0.03 


0.927988 


2.15520 




0.8689 


0.0 


0.928042 


2.15507 




0.8713 


-0.01 


0.92806 


2.15503 


FITO 






0.92819 


2.15473 



TABLE VII: The dominant exponents and the first correction exponents for the relaxation of the 
second moment of magetization with ordered initial states, obtained from the two fitting procedures 
FIT2A and FITO for the case of q = 2. The column for c' M lists the chosen values of the coefficients 
of the second correction terms for fits. 
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